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ABSTRACT 


Many  (perhaps  most)  service  systems,  such  as  repair  and  job  shops, 
computation  centers,  and  transportation  networks,  experience  demand  that 
is  non-stationary  in  time.   This  paper  describes  models  for  situations 
in  which  demands  made  are  by  a  finite  number  of  individuals,  who,  having 
been  served,  do  not  return  until  much  later.   Such  a  transitory  demand 
or  arrival  process  describes  many  phenomena,  among  them  being  commuter 
rush  hours,  and  also  perhaps  the  effect  on  a  population  of  individuals 
their  simultaneous  exposure  to  a  dosage  of  medicine,  a  disease,  or  even 
a  pollutant.   Our  paper  formulates  several  models  for  the  service  of  such 
demands  and  describes  the  manner  in  which  system  state  may  be  approximated 
by  Gaussian  processes,  in  particular  the  Ornstein-Uhlenbeck  and  Wiener 
diffusions . 
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I.   Introduction. 

The  theory  of  service  systems  and  waiting  lines  has  developed 
extensively  over  a  number  of  years,  stimulated  by  a  desire  to  describe 
congestion  and  traffic  in  many  applied  situations.   Almost  without 
exception,  however,  the  mathematical  theory  has  dealt  with  stationary 
models,  in  which  a  steady  stochastic  flow  of  arrivals  approaches  a  set 
of  servicing  facilities;  for  recent  exceptions  see  the  work  of  Newell 
[lCj ,  and  that  of  Milch  [9].  Nevertheless,  many  quite  reasonable  situa- 
tions suggest  the  study  of  models  that  are  non-stationary.   In  this 
paper  we  describe  and  analyze  a  class  of  models  of  a  distinctly  non- 
stationary  nature.  We  call  them  transitory  service  systems. 

Transitory  phenomena  occur  when  members  of  a  finite  population 
of  customers  makes  application  for  service  once,  or  at  most  a  finite 
number  of  times.   Some  examples  follow. 

(a)  A  questionnaire,  letter,  or  newspaper  or  television  advertise- 
ment offering  a  prize  or  job,  is  transmitted.  A  certain  number  of 
individuals  choose  to  respond,  and  do  so  after  various  times.   The 
responses  must  then  be  processed.   Of  interest  is  the  number  of  indi- 
viduals in  various  stages:   e.g.  potential  responders  waiting  to  respond, 


and  responders  undergoing  processing.   One  may  also  wish  to  infer  the 
number  of  potential  responders  in  the  population  on  the  basis  of  the 
responses  during  a  short  initial  time  period. 

(b)  A  message  is  sent  by  radio  to  a  group  of  stations  or  ships. 
The  message  requires  a  reply  or  verification,  an  act  that  takes  place 
after  a  random  time,  and  requires  a  random  time  to  complete.   The 
probability  distribution  of  the  number  of  messages  being  verified  at  a 
particular  moment  after  transmission  may  be  useful,  as  are  other  prop- 
erties of  the  reply  process. 

(c)  Certain  copies  of  a  new  model  of  automobile  or  other  consumer 
(or  industrial  or  military)  product  contain  a  defect  which  is  possibly 
dangerous  and  should  be  modified.   When  the  defect  becomes  evident  a 
recall  or  correction  notice  is  sent  out,  following  which  the  copies 
straggle  in  for  modification.  We  are  interested  in  the  number  of  copies 
that  have  applied  for,  and  are  undergoing,  modification,  at  any  time. 

We  are  also  interested  in  the  length  of  time  to  complete  all  modifica- 
tions. 

(d)  A  ship  puts  out  to  sea  with  all  systems  in  operative  condition. 
Realistically,  however,  each  system  has  a  chance  of  failing  during  a 
voyage.   Let  us  assume,  as  may  hopefully  often  be  reasonable,  that  at 
most  one  failure  per  system  will  occur  during  a  voyage;  it  is  this 
assumption  that  makes  our  process  transitory.  When  the  failure  occurs 

a  repair  process  is  initiated.  We  are  interested  in  total  ship  unavail- 
ability (the  number  of  systems  on  repair)  as  a  function  of  time,  and 
other  related  measures  of  total  system  degradation. 


(e)   It  has  been  suggested  [2]  that  a  road  or  freeway  network  may 
be  considered  to  be  analogous  to  an  infinite  server  queueing  system, 
with  drivers'  times  on  the  freeway  identified  to  be  service  times. 
Clearly,  there  will  be  certain  intervals  (rush-hour  periods)  during 
which  a  large  segment  of  the  population  of  automobiles  will  enter  the 
freeway.   They  will  not  return  until  the  next  rush  hour.   Consequently, 
freeway  occupancy  appears  to  be  a  succession  of  transitory  problems, 
one  for  each  rush-hour  period. 

Other  examples  will  occur  to  the  reader.   We  shall  now  embark 
on  an  analysis  of  the  stochastic  models  suggested.   Our  particular 
emphasis  will  be  upon  approximations  of  a  simple  type. 


II.  "Sufficient  Server"  Models. 

Let  there  be  N  potential  applicants  for  service  in  a  population. 
The  time  of  application  of  customer   i,   measured  from  some  initial 

instant,  is  denoted  by  T..   The  service  time  of  customer   i  is  S.. 

J        x  1 

We  assume  that  at  application  each  customer  may  immediately  begin 
service,  i.e.  there  are  "sufficient  servers."  This  is  analogous  to 
the  infinite  server  model  of  conventional  theory. 
Clearly  applicant   i  is,  at  time  t, 

a)  not  yet  in  service  if  T.  >  t, 

b)  in  service         if  T.  <  t,   but  T.  +  S,  >  t, 

1    '        x    i    ' 

c)  finished  with  service  if  T.  +  S.  <  t. 

1    i 

Consider  the  following  models. 

Model  1.   {T.}   and   {S.}   independently  and  identically  distributed 
(i.i.d.). 
Let  T.   have  the  distribution  F,   and  S   the  distribution 
G;   let  F(t)  =  1  -  F(t) ,   and  G(t)  =  1  -  G(t) .   Then  if  A(t)   denotes 
the  number  of  applicants  who  have  not  yet  applied,  Q(t)   denotes  those 
who  have  applied  and  are  in  service,  and  C(t)   denotes  the  number  of 
those  who  are  finished  at   t  we  can  write  down  directly,  on  the  basis 
of  independence,  the  joint  generating  function 

E[uA(t)  vQ(t)  wC(t)]  =  {uF(t)  +  vG*F(t)  +  wG*F(t)}N.        (2.1) 


where 


G*F(t)  = 


[1  -  G(t-x)]dF(x)  (2.2) 


0 


and 


G*F(t)  = 


G(t-x)dF(x).  (2.3) 


0 


If  we  put  u  =  w  =  1  we  see  that  Q(t),   the  number  in  service  has 
the  generating  function 

E[vQ(t)]  =  {[1  -  G*F(t)]  +  vG*F(t)}N  (2.4) 

of  the  binomial  distribution.   Therefore  the  expectation  and  variance 
of  the  number  in  service  at  time   t   are 


E[Q(t)]  =  N 


t 
f 


[1  -  G(t-x)]dF(x) 


Var[Q(t)]  =  N    [1  -  G(t-x)dF(x){l  -    [1  -  G(t-x) ]dF(x) } . 
0  0 


As  N  becomes  large  one  can  approximate  the  distribution  of  Q(t)   by 
the  normal  with  mean  and  variance  matching  (2.5). 

The  conventional  infinite  server  problem,  to  which  this  model 
compares,  yields  no  such  simple  time-dependent  result.   Recall  that  if 
arrivals  are  stationary  Poisson  then  the  stationary  distribution  is 
also  Poisson;  see  Riordan  [12].   Of  course,  if  population  size,   N, 
is  large  then  Q(t)   may  turn  out  to  be  approximately  Poisson  at  a 
fixed  time   t.   Consider  the  following  extensions. 


Model  2.   {T  }   and   {S  }   i.i.d.;   N   is  Poisson. 

If  population  size  N   is  a  random  variable  with  the  Poisson 
distribution  of  mean  a,   then  clearly  Q(t)   is  also  Poisson  with 
mean  aG*F(t) ,   as  follows  from  the  generating  function: 


oo 

E[vQ(t)]  =      I      {I   -   G*F(t)  +  vG*F(t)}n 
n=0 


=  exp{aG*F(t)(l-v)}. 


n 
-a  a 

— r 
n! 


(2.6) 


Moreover,  A(t) ,   Q(t),   and  C(t)   are  jointly  and  independently 
Poisson  in  this  model  as  is  evident  from  the  fact  that  the  generating 
function  factors: 

E[uA(t)vQ(t)wC(t)]  =  exp{aF(t)(l-u)}exp{aG*F(t)(l-v)} 

(2.7) 
exp{aG*F(t)(l-w)}. 

It  cannot,  however,  be  concluded  that   {Q(t),t^0}   is  a  Poisson  process 


Model  3.   {T.}   are  i.i.d.  Gamma;   {S.}   are  i.i.d.;   the  population 
size  N  becomes  large. 
Suppose  the  arrival  distribution,  that  of  T.,   is  Gamma,  with 
density 

krl 


f(t)  =  e  Xt^|}X  (2.8) 


(k) 

for   t  ^  0  and  k  >  1.   Then  we  can  express  the  generating  function 
of   Q(t),   (2.4),  as 


E[vQ(t)]  -  [1  -  (1-v) 


[l-Gtt^le-^l^g-^dxf     (2.9) 


Now  scale  time  by  letting  X  N  =  6,   a  constant,  and  allow  N  -»■  ». 


Then 

t 


E[vQ(t)]  -+  exp{-(l-v)6 


k-1 
[1  -  G(t-x)]  Xr(k)dX  }       (2.10) 


0 


t  k-1 

[1  -  G(t-x)]      "  .  Notice  that  if  k  =  1  then  in 
0  Iu; 


The  continuity  theorem  then  implies  that  under  the  stated 
conditions  the  distribution  of  Q(t)  approaches  the  Poisson  with 
paramters  6 

the  limit  we  have  the  familiar  "infinite  server"  solution  of  the  classi- 
cal theory:   if  arrivals  are  a  stationary  Poisson  process — to  which 

the  process  {A(t),t^0}   converges  as  N  ->  °° — then  {Q(t),t^0}   is  a 

•t 
non-stationary  Poisson  process  with  E[Q(t)]  =  9     [1-G(x)]dx  ->  6E[S] 

as   t  -*■  oo.   For  other  k  values  the  process  behaves  in  a  comparable 

fashion.   It  may  be  contended  that  the  type  of  large-N  result  just 

derived  actually  justifies  the  use  of  the  classical  queueing  models 

that  assume  Poisson  arrivals  with  stationary  independent  increments, 

at  least  early  in  the  arrival  process.   Of  course,  an  approximation 

based  upon  NX  =  9   cannot  be  accurate  late  in  the  process,  e.g.  after 

In   2 
a  time  t.  ,_  =  — - — ,  when  on  the  average  one-half  of  the  calling 

population  has  arrived.   Therefore,  in  the  following  section  we  develop 

a  procedure  for  approximating  our  processes  A(t) ,  Q(t),   and  C(t)   by 

Ornstein-Uhlenbeck  processes;  these  approximations  can  be  expected  to 

be  reasonably  accurate  for  moderate  N,   and   t-values  of  practical 

interest. 


III.   Diffusion  Approximation  for  Sufficient  Servers. 

Since  A(t) ,   the  number  of  arrivals  in  time  t,   and  Q(t), 
the  number  in  service,  are  approximately  normally  distributed  for  large 
N,  we  are  led  to  consider  the  use  of  a  diffusion  approximation  for 
these  processes;  see  Feller  [4],  and  Cox  and  Miller  [3].   We  remind  the 
reader  that  a  diffusion  is  a  Gaussian  process  with  continuous  sample 
paths.   We  next  proceed  to  develop  approximating  diffusions  for  certain 
transitory  systems,  and  then  to  make  use  of  such  approximations  to 
obtain  properties  of  the  service  systems.   Our  approach  here  is  in  the 
spirit  of  McNeil   [8];  other  authors  have  utilized  similar  procedures. 
For  example,  see  Anderson  and  Darling  [1],  and  Whittle  [14]. 

III-l.  Approximating  the  Arrival  Process. 

Let  \j(t)   denote  the  number  of  arrivals  by  t  out  of  a  popula- 
tion of  N.   Each  arrival  time  is  independently  distributed  with  d.f . 
F(t),   density  function  f (t) ,   and  hazard  rate  X(t)  =  f (t) [l-F(t) ]~  . 

Now  clearly  — ►  F(t)   in  probability  as  N  •*   °°  for  every 

fixed  t.   In  fact,  the  convergence  is  uniform  with  probability  one  by 
the  Glivenko-Cantelli  theorem;  see  Tucker  [13],  p.  127.  Hence,  we  are 
led  first  to  approximate  \,(t)   by  a  deterministic  component,  x(t) , 
for  large  N: 

^(t)  ~Nx(t)  (3.1) 

Noting  that  if  \r(t)   arrivals  have  occurred  by  t,   and  thus  that  the 
probability  of  another  arrival  in  (t,t+dt)   is  essentially 
[N  -  A^(t)]X(t)dt  we  can  derive 

^=  X(t)[l  -  x(t)],  (3.2) 

the  solution  of  which  is  x(t)  =  F(t). 


To  include  a  stochastic  element,  or  noise,  study 

Aft)  -  Nx(t) 

S  (t)  =  -^ .  (3.3) 

N  v¥ 

Our  first  approach  is  to  derive  the  transformed  version  of  the  forward 

differential  equation  satisfied  by  the  density  function  of   S  ,   pass 

to  the  limit  as  N  ->  °°,   and  thus  show  that  if   S   converges  to  a 

limiting  process  the  latter  must  be  non-stationary  Ornstein-Uhlenbeck; 

see  Cox  and  Miller  [3] . 

If 

i6AN(t)  r- 

♦N(6,t)  =  E[e   W   ]     i  -  /=T,  (3.4) 

is  the  characteristic  function  (ch.f.)  of  \,(t)  we  have 

i9Aft+dt)       ie\.(t) 
E[e   N      ]  =  E[e   W   {1  -  [N-A^t)  ]X(t)dt}] 

i6(A(t)+l) 
+  E[e  {N-AN(t)}X(t)dt]  +  o(dt),    (3.5) 

and  upon  collecting  terms  and  letting  dt  -*■  0   there  emerges  the 
differential  equation 

—f^   =  (eie-l){N^(e,t)  +  i  -1f>X(t)  (3.6) 

Now  to  derive  an  equation  for  the  ch.f.  of   S  (t) , 

ies  (t) 

*N(t)  =  E[e   N    ],  (3.7) 

we  note  that 


(i)  ^(e.t)  -  e-iex(t)vStlI(e/^.t), 


(11)  ^e/^.t)  -  .««<t)^!!H^!i+  ux-w^r  .»«<«>*  #,,(6.0. 
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(iii)  ^f(e/^,t)  -  .«*<«■*  ^e.t)  +  ix(t)^e19x(t>^^(6,t).   (3.8) 

Next  change  6   to  6'  =  Q/vft     in  (3.5),  utilize  (3.7),  expand  the 

exponentials   (e      -1= on"  +  Q(  o/9))   and  collect  terms.   The 

essential  outcome  is 

8*N       jc-  62 

-^+  ^Ni6v¥{x'(t)  -  X(t)[l-x(t)]}  =  -|-X(t)[l-x(t)]^N 

3*N      1 
-  6X(t)  -r£  +  o(-±-).        (3.9) 

In  order  for  the  equation  to  be  satisfied  by  a  limiting  ch.f.,   if(6,t), 
as  N  -*■   °°  the  term  in  brackets  must  vanish;  solution  of  the  resulting 
differential  equation  once  again  presents  the  deterministic  equation 
of  (3.2).  There  remains  the  equation 

|i  .  .  ,x(t)  |i  -  £  f(e),  O.io) 

which  is  the  transformed  version  of  the  forward  equation  for  a  non- 
stationary  Ornstein-Uhlenbeck  (O.U.)  process;  see  Cox  and  Miller  [3], 
pp.  218-219.  The  infinitesimal  mean  and  variance  are  respectively 

3(x,t)  =  -xX(t)   and  a(x,t)  =  f(t).   Since  the  O.U.  process  is  Gaussian 

82 
we  may  differentiate  the  Gaussian  ch.f.   exp[i6u(t)  — ~—  c2(t)]   and 

equate  coefficients  in  (3.9)  to  obtain  the  mean  and  variance.  Not 

surprisingly, 

E[S(t)]  =  0,     and  Var [S (t) ]  =  F(t) [l-F(t) ]        (3.11) 

Hence  we  are  led  to  approximate  \r(t)   by 

^(t)  «NF(t)  +  /N  S(t)  (3.12) 

where  S(t)   is  an  O.U.  process  with  parameters  given  by  (3.11). 


11 


An  alternative  derivation  of  the  O.U.  parameters  is  easy,  once 
one  settles  upon  a  diffusion  approximation.   The  stochastic  differential 
equation  for  A^   is  given  approximately  by 

dAN(t)  =  X(t)[N-AN(t)]dt  +  Z(t)A(t)[N-AN(t)]dt,       (3.13) 

Z(t)   being  a  "purely  random  process,"  or  the  "derivative"  of  a  Wiener 
process.   The  first  right-hand  term  of  (3.12)  is  the  deterministic  com- 
ponent or  drift,  and  the  second  is  the  supplementary  randomness  or  noise, 
the  scale  of  which  follows  from  considering  possible  events  in  a  dt 
interval.   Now  apply  normalization  (3.3)  to  obtain 

dSN  +  v¥[x'(t)  -  X(t){l-x(t)}]  =  -  X(t)SN(t)dt 


+  Z(t)/X(t 


i-x(t)  -  ^ 


dt,  (3.14) 


and  once  again  we  are  compelled  to  set  x'(t)  =  X(t)[l-x(t)]   and  to 
realize  that  the  limiting  noise  satisfies 

dS(t)  =  -  X(t)S(t)  +  Z(t)/f(t)dt,  (3.15) 

the  stochastic  differential  equation  for  the  O.U.  process  with  drift 
and  variance  parameters  already  derived. 

111-2.   Passage  Probabilities  for  the  Arrival  Process. 

In  this  section  we  show  how  the  limiting  O.U.  process  furnishes 
an  approximation  to  the  probability  that  \,(t)  £  8(t)   for  all   t,   8(t) 
being  a  suitably  selected  boundary.   Our  asymptotic  result  is  closely 
related  to  the  one-sided  Kolmogorov-Smirnov  test,  and  may  be  used  for 


12 


the  same  purpose,  i.e.  as  a  test  of  the  statistical  hypothesis  that 
the  arrival  distribution  F(t)   generates  arrivals  more  rapidly  than 
does  F  (t)  ,   i.e.  that  F(t)  2>  F  (t)   for  all   t.   The  methodology 
employed  may  also  be  used  to  estimate  the  maximum  queue  size  during  a 
single  rush  period,  but  further  approximations  are  required  and  this 
problem  will  be  attacked  later. 

It  is  shown  in  [3]  that  an  O.U.  process,   S(t)  may  be  expressed 
in  terms  of  a  Wiener  process,   X(t) : 

S(t)  =  a(t)X{x(t)}  (3.16) 

where  x(t)   represents  a  time  scale  transformation,  and  a(t)   is  a 
real  dif ferentiable  function.   In  order  to  represent  our  particular 
process  we  have 

|^-=-A(t),   and   a2(t)x'(t)  =  f(t).         (3.17) 

Solving,  we  find  that 

a(t)  =  1  -  F(t),   and   x(t)  =  ±  l^\t)    •         (3.18) 

Hence 

X(t)   being  a  Wiener  process  with  zero  drift  and  infinitesimal  variance 
t.   It  follows  that  if 

X(t)  £  8(t),   for  {t  :  0  £  x  £  F(T) [l-F(T) ]_1} 


then         S(t) 

±  _l[t)   *   B(x(t))   for   {t:0^t^T}. 
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We  consider  certain  special  cases  results  for  which  are  at  hand, 
Linear  Boundaries.   Let   8(t)  =  a  +  3t,  Vt  ;>  0;   a, 3  >  0.   Then 
it  is  known  that 


P{X(t)  £  a  +  8t,  Vt  :>  0}  =  1  -  e 


-2aB 


(3.20) 


Hence 


S(t)  £  a[l-F(t)]  +  3F(t),   Vt  :>  0 


_2fY  R 
with  probability  1  -  e    ,   and  furthermore 


Aft)  <;  NF(t)  +  /N{a[l-F(t)]  +  3F(t)},    Vt  ^  0       (3.21) 

_2(Y  ft 

with  approximate  probability  1  -  e 

Linear  Boundaries,  Finite  Stopping.  Let  B(t)   be  as  above;  then  if  X 

is  a  Wiener  process  and  t  >  0  is  fixed, 

P{X(t')  <  a  +  At',  V  t1  <  t,  X(t)  si  x  <  a}  = 


x-gr-i     2a B^  (-x-ftT-2a 
e 


/t 


(3.22) 


where  <J>   is  the  standard  normal  distribution;  see  Cox  and  Miller  [3], 
pp.  220-221.  Now  use  (3.19): 


P{S(t)  s;  a[l-F(t)]  +  3F(t),  V  t  £  T,  S(T)  <:  y}  = 


y(l+T)-0T 


where  t  = 


F(T) 


l-F(T) 


T  =  F 


-1 


e     $ 


y(l+T)-3T-2a 


/F 


(3.23) 
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Thus  we  can  represent  the  probability  that  the  noise  component  at  time 
T  is  in  any  Borel  set  of   (-«,  a[l-F(T)]  +  6F(T))   for  any  finite 
T.  This  result  can  be  carried  over  to  an  approximation  to  the  distribu- 
tion of  A^(t) ,   subject  to  boundary  avoidance  prior  to  T. 

Similar  problems  involving  two  boundaries  are  also  solvable. 
In  addition,  the  distribution  of  the  time  required  to  cross  a  boundary 
of  the  form  a[l-F(t)]  +  3F(t)   can  be  derived,  although  moments  are 
not  usually  expressible  in  a  simple  closed  form. 

Application  of  the  above  material  to  a  problem  in  sequential 
testing  is  described  in  the  thesis  by  Gwinn  [5],  and  an  expanded  treat- 
ment is  in  preparation. 

III-2.   Diffusion  Approximation  of   (A^(t) ,Q  (t)) . 

In  order  to  establish  a  diffusion  approximation  for  QN(t) , 

the  number  undergoing  service,  we  must  treat  the  bivariate  process 

(A^(t),Q  (t)),   focusing  on  the  stochastic  elements 

Aft)  -  Nx(t) 
S„(t)  =  -2- 


N        >¥ 


and 


Q  (t)  -  Ny(t) 

TN(t)  =  -3 (3.24) 

N  vfi" 

We  first  proceed  as  was  done  earlier.  Letting  X(t)  =  f(t)[l-F(t)] 

denote  an  individual's  arrival  rate,  and  assuming  exponential  service 

times,  we  write 

P{AN(t)  =  n,QN(t)  =  jlAjjCt)  -  n,QN(t)  =  j }  =  l-{X(t)  [N-n]+uj}dt+o(dt) 
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P{AN(t)  =  n+l,QN(t)  =  j+llA^t)  =  n,QN(t)  =  j}  =  X(t)  [N-n]dt+o(dt) 

p<AN<t)  =  n»QN(t)  =  J-llA^t)  =  n,QN(t)  -  j}  -  yjdt+o(dt).     (3.25) 

other  probabilities  being  negligible. 

Next  form  the  characteristic  function  expression 

ie  A^(t+dt)+ie  q  (t+dt)     ie  iL(t)+ie  q  (t) 

E[e  ]  =  E[e   X  N      2N    {1-X  (t)  [N-A^t)  ]  }x 

ie  +ie  -ie 

{l-e  X   Z}  -  uQN(t)U-e   Z}  +  o(dt)    (3.26) 

Then  substitute  from  (3.24),  divide  by  dt  ->  0,   and  substitute  as  in 
(3.8),  letting  N  ->  «>.   in  order  for  a  limiting  ch.f.  \\>     to  exist  for 
(S  ,T  )   it  turns  out  that 

x*(t)  =  X(t)[l-X(t)]  =  f(t) 
and  (3.27) 

y'(t)  =  X(t)[l-x(t)]  -  uy(t) 

and  solution  gives  the  deterministic  components  of  the  previous  section. 
Finally,   ip   satisfies 

It = "  i[(ei+e2)2  f(t)  +  ^y^*  -  (Q1+Q2n(t)  w~  * "  V  a!-  *     (3-28) 

This  partial  differential  equation  is  the  Fourier  transformed  version 
of  a  bivariate  non-stationary  Ornstein-Uhlenbeck  process. 

A  stochastic  differential  equation  approach  runs  as  follows. 
Approximately, 


dQN(t)  =  {X(t)[N-AN(t)]  -  uQN(t)}dt  +  Z(t)/u(t)[N-AN(t)]dt+  QN<t)dt,     (3.29) 
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so  when  (3.24)  is  applied  and  N  -+  °°  we  find  that  the  deterministic 
part,   y,   must  satisfy 


y1  =  f(t)  -  uy 
or  (3.30) 


y(t)  = 


-u(t-s)  ,M, 
e       f(s;ds 


0 
in  accordance  with  earlier  results,  while  the  limiting  noise  is  described 

by 


dT  =  -[X(t)S(t)  +  uT(t)]dt  +  Z(t)/[f(t)+uy(t)]dt      (3.31) 

Thus  T   is  conditionally  Ornstein-Uhlenbeck,  given  S(t).   Also,  the 
purely  random  parts  of  dT  and  dS  are  correlated,  as  is  clear  from 
the  fact  that  an  additional  arrival  lengthens  the  queue;  the  correlation 
term  is  f(t),  and  we  can  write 

dS(t)  =  -X(t)S(t)dt  +  Z  (t)/f(t)dt 

dT(t)  =  -[uT(t)+A(t)S(t)]dt  +  Zx(t)/f(t)dt  +  Z2(t)/uy(t)dt    (3.32) 

where  Z.   and  Z   represent  independent  purely  random  components,  i.e. 
Wiener  process  derivatives. 

I1I-3.  Moments. 

Since  {S(t),T(t)}   is  a  bivariate  diffusion  with  zero  mean  we 
know  that  ^  is  of  the  Gaussian  form 

^(91,e2,t)  =  exp{-  j  6/  Z_0_}  (3.33) 

where  9/  =  (0- ,0_) ,   and  2^  is  the  covariance  matrix.  Hence,  we  need 
merely  equate  coefficients  of  0  9    (i,j  =  1,2)   in  order  to  derive 
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the  covariance  matrix  elements.   Not  surprisingly,  these  agree  with 
earlier  exact  results: 

Var[S(t)]  =  x(t)[l-x(t)];  Var[T(t)]  =  y(t)[l-y(t)] 
E  F(t)[l-F(t)] 


Cov[S(t),T(t)]  =  [l-F(t)] 


-u(t-s)  ,,  NJ 
e        f(s)ds. 


(3.34) 


0 
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IV.   Single-Server  Models. 

Many  of  the  transitory  situations  that  are  encountered  in 
practice  actually  involve  competition  for  facilities;  e.g.  for  repairmen, 
clerks,  freeway  space,  or  computers.  Therefore  in  this  section  we 
describe  models  for  a  single-server  system  confronted  by  transitory 
demand . 

Once  again,  a  bivariate  stochastic  process  is  required  to  describe 
our  model.   If  N   individuals  wait  to  make  application  for  service, 
then  the  number  of  new  arrivals  in  (t,t+dt)   depends  upon  the  random 
number  that  has  already  appeared  by  time  t.   The  waiting  line  develop- 
ment is,  of  course,  influenced  by  the  service  rate  as  well.   It  is  the 
coordinated  effects  of  these  processes  that  we  represent,  and  attempt 
to  make  comprehensible  by  means  of  approximations . 

Model  1.  Markovian  Arrivals  and  Exponential  Service. 

Let  A(t)   be  the  number  of  arrivals  that  have  occurred  by  t, 
and  Q(t)   be  the  number  waiting  and  in  service  at  t  at  the  service 
facility.   Then  let 

P.  (t)  =  P{Q(t)  =  j,  A(t)  =  n  |  A(0)  =  0,  Q(0)  =  0},     (4.1) 

for  0  £  j  s:  n  £  N.   Given  A(t)  =  n,   let  A  dt  be  the  probability 
of  exactly  one  arrival  in  (t,t+dt),   and  let  udt  be  the  probability 
of  exactly  one  departure  in  the  same  interval.  The  usual  independence 
assumptions,  cf.  Feller  [4],  are  made.   Then  we  can  write  down  the 
forward  differential  equations  for  {P.  (t)}: 
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poo(t)  =  "  xo  poo(t) 


Pin(t)  =  "(Xn+M)Pin(t)  +  ^+1  n(t>  +  X   1  P-  1     1  <*>  • 

jn        n    jn       j+l,n       n-1  j-l,n-l 


j  =  1,2,. . .,n-l; 


P0n(t)  =  "  Xn  P0n(t)  +  *  Pln(t) 


P'  (t)  =  -(X  +u)P   (t)  +  X  .  P  -    n(t),  (A. 2) 

nn        n    nn       n-1  n-1, n-1 

where  initial  conditions,  e.g.   P~n(0)  =  1,   and  P,  (0)  =  0  otherwise 

UU  jn 

must  be  specified.   The  arguments  leading  to  (4.2)  are  standard  and 
will  be  omitted.   Of  course,  one  can  also  allow  arrival  and  departure 
rates  to  depend  upon  t. 

For  any  explicitly  specified   X    the  solution  to  (4.2)  can 
be  computed  numerically,  for  (4.2)  is  only  a  system  of  finitely  many 
differential  equations.   If  desired,  one  can  Laplace  transform  through- 
out, but  there  is  no  simple  explicit  formula  even  for  the  transform. 
See  the  paper  by  Perlas,  [11]  ,  for  further  details. 

Model  1-A.   Independent  Exponential  Arrivals. 

A  model  of  some  interest  assumes  that  X     =  X (N-n) ,   or, 

n 

equivalently,  that  each  individual's  arrival  time  is  independent  with 
distribution  1  -  e  .  Then  a  deterministic  approximation,  that  one 
would  be  tempted  to  use  for  N  large,  is  suggested: 

dx(t)  =  X[N-x(t)]dt  (4.3) 

and 

y(t+dt)  =  max{y_(t)  +  X[N-x(t)]dt  -  ydt.O},         (4.4) 
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where  x     and  y_  represent  deterministic  approximations  to  A  and  Q, 
Solutions,  subject  to  A(0)  =  Q(0)  =  0,   are 

x(t)  =  N[l-e"Xt] 


and 


y_(t)  =  max{N(l-e~Xt)  -  ut,0}.  (4.4) 


More  will  be  said  about  these  approximations  in  the  next  section. 

Numerical  solutions  of  the  above  equations  have  been  carried 
out  to  compare  the  total  waiting  times  (i)  under  the  full  stochastic 
model,  (4.1),  and  (ii)  under  the  deterministic  approximation  (4.2),  (4.3). 
These  are  summarized  in  Tables  1  and  2,  which  show  that  the  deterministic 
approximation  improves  greatly  as  N  ->  °°.  Not  surprisingly,  the  deter- 
ministic approximation  falls  below  the  stochastic  expectation. 

IV-1.   The  Probability  of  No  Wait. 

In  [6],  Kabak  describes  a  transitory  problem  in  which  atten- 
tion focuses  upon  the  probability  that  no  arrival  must  wait.  Other 
applications  for  this  model  may  well  exist,  e.g.  in  military  situations. 
Let  us  suppose  that  calls  occur  in  accordance  with  the  pure  birth  process 

of  Model  1,  and  let  successive  service  times   {S  ,   n  =  1,2,...N}   be 

n 

independent  but  arbitrarily  distributed  in  accordance  with  F  (x) .   If 

T,.s      represents  the  time  until  the  i —  arrival  occurs,  then  the 
(i) 

probability  of  no  wait  is  seen  to  be  that  of  the  following  event 

T(2)>S1,T(3)"T(2)>S2»T(4)"T(3)>S3,"*,T(N)"T(N-1)>SN-1    (4*5) 

Since  times  between  successive  arrivals  are  independently  exponential 

in  this  Markovian  model  it  is  immediate  that  the  probability  of  no  wait  is 


N-l 

p  =  n 
j-i 


00 

t 


0 
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-X  t  N-l 

e  J   dF  (t)  =  n   F  (X  )  (4.6) 

J      j=1  J   J 


where  F   is  the  Laplace-Stieltjes  transform  of   F. 
In  the  case  of  Model  1-A  this  specializes  to 


p 

— 

V 

y 

►  •  • 

y 

y  + 

Xl 

y  +  x2 

y 

+  Vi 

N- 

-i 

£ 

+  N- 

-1)  (*  +  >. 

-2) 

• 

..  (*+ 

i) 

(5 

N 

:)   r® 

(N- 

■D! 

y 

(N-l)X 

as      N  -> 

(A. 7) 


(4.7) 


according  to  Euler's  product  definition  of  the  gamma  function;  see  Knopp 
[7].   Of  course,  if  service  times  are  constant  rather  than  exponential 
a  simple  exact  solution  is  possible: 

N-l 
-;  I  X(N-j)     _  X  N(N-l) 

P=e    J=1       =e"^   2 


In  the  event  that  F   is  a  one-sided  stable  distribution  of  order  a, 

0  <  a  <  1,  we  have  essentially 

N"1  i, 

-Id  l   Cn-J)  -(h    N —         (4  9) 

y   j=l  V   1+a             V   ; 

P  =  e      J  ~  e 


Thus  even  though  the  stable  laws  have  right  tails  long  enough  to  engender 
infinite  first  moments,  the  mass  near  zero  governs  the  above  probability. 
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IV-2.  The  Duration  of  the  Process. 

If  we  let  t  denote  the  time  to  complete  the  process  of  arrival, 
queueing,  and  service,  then  it  is  clear  from  (4.1)  that 

P{t  ac  t}  -  PQN(t).  (4.10) 

Unfortunately,  this  distribution  has  no  neat  closed  form,  although 
numerical  properties  can  be  derived  by  solution  of  (4.2) . 

It  is  worth  noting  that  for  large  N  and  under  the  assumption 

of  Model  I-A  (4.10)  is  very  close  to  a  gamma  distribution  with  mean 

N  N 

—  and  variance  — z".   The  explanation  is  that  for  large  N  a  queue 

quickly  develops  behind  the  server,  and  does  not  vanish  until  the  last 

service  occurs.   The  server  effectively  continues  to  be  active  through 

the  service  of  all  N  arrivals. 
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V.   Diffusion  Approximation  for  Single-Server  Models. 

The  setup  to  be  analyzed  now  is  as  follows .   N   individuals 

independently  choose  their  arrival  times  from  a  distribution  F(x) , 

the  latter  possessing  a  smooth  density  f(x).   In  the  order  of  their 

arrival  these  are  served  at  a  single  servicing  facility;  service  times 

are  exponential  at  rate  Nu .   Let  \,(t)   denote  the  number  of  arrivals 

that  have  occurred  on  or  before   t,   QvjCO   i-s  tne  number  awaiting  or 

•t 


undergoing  service  at   t,   and  D N(t)  = 


Q  (s)ds   is  the  total  delay 
0 


accumulated  up  to  time   t.   We  shall  analyze  the  birth  and  death  process 
{A^(t) ,Q  (t) } ;   if  F   is  exponential  this  is  just  the  process  discussed 
in  the  last  section.   A  useful  approximate  approach  is  to  derive  a 
diffusion  approximation,  valid  for  large  N.   In  the  present  approxima- 
tion,  F(t)  >  ut   for  a  substantial  range  of   t-values  (heavy  traffic) 
is  also  a  requirement. 

V-l.   Deterministic  Approximation  of   {A^(t) ,Q  (t) ,D  (t) } . 

Begin  by  isolating  the  deterministic  component  of  the  process 

Q  (t) ,   since  that  of  \,(t)   has  already  been  treated.   Notice  that  if 

QN(t) 

— >-  y(t)   as  N  ->-  a>  and  it  is  assumed  that  y  has  a  derivative, 

then 

QN(t+h)  -  QN(t)  =  X(t)[N-AN(t)]dt  -  Nydt  +  o(dt)        (5.1) 

if  Q  >  0,   while 
N 

QN(t+h)  -  QN(t)  =  X(t)[N-AN(t)]dt  +  o(dt)  (5.2) 

if  Q  =0.   Dividing  by  dt,   and  by  N,   and  taking  the  limit  as   dt  ■>  0 
and  N  ->  oo  yields  the  differential  equation 
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y'(t)  =  X(t)[l-x(t)]  -  uh[y(t)]  (5.3) 

where  y(0)  =  0  and 

{1,     x  >  0 
(5.4) 
0,     x  £  0; 

the  solution  is 

y(t)  =  max{F(t)  -  yt,0}.  (5.5) 

It  also  follows  that   z(t),   the  deterministic  component  of 
D(t),   is 


z(t)  = 


y(s)ds  (5.6) 


0 


V-2.  Diffusion  Approximations. 

Next  define  the  supplementary  randomness  or  noise  components 
of  our  three  processes: 

Aft)  -  Nx(t) 

s»(t)  =      \^ ' 

Q  (t)  -  Ny(t) 

T„(t)  =  -H (5.7) 

N  /S 

D  (t)  -  Nz(t) 
UM(t)  -  — 

N        v¥ 

Begin  by  writing  down  the  infinitesimal  transition  probabilities: 
P{AN(t+dt)  =  A^t)  ,  QN(t+dt)  =  QN(t),DN(t+dt)  =  DN(t)+dtQN(t)| 
AN(t),QN(t),DN(t)}  =  l-X(t)[N-AN(t)]dt-yNh[QN(t)]dt+o(dt), 
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P{AN(t+dt)    =  AN(t)+l,QN(t+dt)    =   QN(t)+l,DN(t+h)    =   DN(t)+dt[QN(t)+|    | 

AN(t),QN(t),DN(t)}    =   X(t)[N-AN(t)]dt+o(dt) 

(5.8) 
PCAjjCt+dt)    =  AN(t),QN(t+dt)    =    [QN(t)-l]h[QN(t)],DN(t+dt) 

=  DN(t)+dt[QN(t)-|]     |    AN(t),QN(t),DN(t)}      =  UNh[QN(t)]dt+o(dt) 

where  once  again  h(«)   represents  the  unit  step  at  the  origin.   Now 
put  for  the  joint  ch.f. 

ie  s  (t)+ie  t  (t)+ie  u  (t) 

*N(61,82,e3;t)  =  E[e  L  N      l   N  ]       (5.9) 

In  order  to  derive  a  partial  differential  equation  for  \\>        follow  the 

conditioning  argument  used  earlier  on  \,(t)   in  section  III.   Also, 

assume  throughout  that  h[Q  (t)]  =  1,   so  our  approximate  solutions  will 

not  apply  when  Q   is  small  (we  must  avoid  the  very  end  of  the  process, 

and  for  some  distributions,   F,   the  very  beginning).   Omitting  the 

tedious  algebra  we  state  that  if  a  limiting  noise  distribution  exists 

its  ch.f.  \\>  =   lim  \p    ,   must  satisfy  the  differential  equation 
N-*» 


||  =  -  |[f(t)(e1+e2)2  +  \iQ*U  -  x(t)(e1+e2)  ^|-*+  e3  g|- *   (5, 


10) 


This  equation  is  the  Fourier-transformed  forward  equation  of  a  trivariate 
non-stationary  diffusion.   As  we  already  know,   S(t)   is  non-stationary 
O.U.   Conditionally  upon  S(t),   T(t)   is  a  non-stationary  Wiener  process, 
This  is  independently  established  by  deriving  the  stochastic  differential 
equation: 
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dQN  =  U(t)[N-Nx(t)-vfiSN(t)  -  Ny}dt 


+  Z(t)/{X(t)[N-Nx(t)-^SN(t)+Nu  dt      (5.11) 


=  /N  dT.T  +  Ndy 

N 

After  passage  to  the  limit  as  N  -*■  °°  and  excision  of  the  deterministic 
component  the  expression 

dT  =  -X(t)S(t)dt  +  Z'(t)/{f(t)+y}dt  (5.12) 

results.   We  note  that  (5.10)  may  be  derived  directly  from  (3.15)  and 
(5.12),  the  stochastic  differential  equations  for  the  arrival  and  queue 
noise  components.   Forgetting  about  the  U  component,  write 

^(e1,e2,t+dt)  -  i^(e1,e2,t)  = 

ie1{s(t)+ds>+ie2{T(t)+dt}       ie1s(t)+ie2T(t) 

E[e  -  e  ] 


ie^s+ie  dt     ie1s(t)+ie9T(t) 

=  E[{e  L  -   l}{e  }]. 


Now  make  use  of  (3.15)  and  (5.12),  and  the  additional  fact  that  the 
purely  random  components  are  correlated: 


^(e1,e2,t+dt)-i|)(e1,e2,t)  =  E[{e 


-i61X(t)S(t)dt-ie2X(t)S(t)dt 


^  {f(t)(e1+e9)2+ye9}dt  ie.s(t)+ie9T(t) 

1  \ii     _1}x  e    i  *■       j      (5#14) 
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Next  expand  exponents  to  order  dt,   divide  by  dt,   and  let  dt  -►  0; 
after  a  little  rearrangement  (5.14)  turns  into  (5.10),  with   6=0 
and  hence  the  last  term  removed.   Note  that  the  purely  random  parts  of 
S   and  T  are  correlated  with  covariance  equal  to   f(t)dt.   Thus  one 
might  write  the  equations  as  follows: 


dS  =  -A(t)S(t)dt  +  Z  (t)/f  (t)dt 

dT  =  -X(t)S(t)dt  +  Z  (t)/f(t)dt  +  Z  (t)/pdt         (5.15) 

where  Z   and  Z   are  independent  and  purely  random.   It  now  follows — 
as  it  does  also  by  setting   9.  =  -Q    t      and   6„  =  0   in  (5.10) —  that 
the  noise  process  T(t)  -  S(t)   is  a  Wiener  process  with  zero  drift  and 
infinitesimal  variance  u. 

In  summary,  our  assumptions  (which  are  akin  to  those  of  heavy 
traffic  theory  in  the  sense  that  boundary  effects  are  essentially  ignored) 
imply  that  arrivals,  number  in  system,  and  total  delay  are  jointly  non- 
stationary  Gaussian.   The  actual  covariance  function  will  be  derived 
subsequently,  and  some  further  approximations  suggested.   From  (5.15) 
it  is  easily  seen  that  a  simulation  of  our  approximation  can  be  carried 
out  by  discretizing  time:   if   t  =  nA,   n  =  1,2,...,   then 


S((n+1)A)  =  S(nA)  -  X(nA)S(nA)A  +  Z  (n)/f(nA)A 


T((n+1)A)  =  T(nA)  -  X(nA)S(nA)A  +  Z  (n)/f(nA)A  +  Z2(n)/yA     (5.16) 

where  here   {Z  (n)  }   and   {Z  (n)  }   are  mutually  independent  sequences 
of  independent  N(0,1)   random  variables,  and  S(0)  =  T(0)  =  0.   Notice 
too  that  the  derivation  of  (5.12)  may  be  carried  out  if  the  service  rate, 
u,   is  a  sufficiently  smooth  function  of   t. 
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V-3.   Moments. 

The  diffusion  process  characterized  by  (5.10)  is  zero  mean 
trivariate  Gaussian,  so 

^(e1,e2,e3,t)  =  exP{-  j  e.'z_  e_}.         (5.17) 


where  9_'  =  [0.,  ,0«t0~] ,   0_  its  transpose,  and  Z_ 


the  covariance 


matrix.   Now  if  we  differentiate  through  with  respect  to   t  and  9  , 
in  accordance  with  (5.10)  and  equate  coefficients  of   9.6.,  we  obtain 
the  following  moments 

Var[S]  =  F(t)[l-F(t)],   Var[T]  =  F(t)[l-F(t)]  +  ut 

Cov[S,T]  =  F(t)[l-F(t)].    (5.18) 

Even  more  easily,  these  moments  can  be  obtained  from  first  principles, 
making  use  of  the  fact  that  Var[T-S]  =  ut  as  implied  by  (5.15).  The 
fact  that  the  boundary  at  zero  has  been  neglected,  and  hence  that  our 
model  is  inadequate  late  in  the  process,  becomes  evident  in  the  expres- 
sion for  Var[T],  which  should  apparently  approach  zero  with  large  t, 
but  does  not.  The  model  should  however  be  accurate  near  the  peak  of 
the  process. 
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VI .   Conclusion. 

The  purpose  of  this  paper  is  to  formulate  and  explore  a  new  class 
of  non-stationary  service  models.   We  have  shown  that  under  certain 
conditions  our  models  may  be  described  by  Gaussian  diffusion  processes, 
and  have  hinted  at  the  manner  in  which  known  results  for  such  diffusion 
may  be  utilized  to  study  model  properties;  for  examples  we  exhibit 
boundaries  below  which  the  process  moves  with  prescribed  probability, 
and  derive  equations  for  the  total  time  spent  waiting  by  all  arrivals. 
Although  all  technical  details  of  approximation  accuracy  are  not  yet 
well  understood,  further  theoretical  and  numerical  explorations  are 
under  way.   In  any  case,  the  diffusion  approximations  suggested  appear 
to  make  an  offer  that  the  modeler  of  stochastic  phenomena  can't  refuse. 

A  brief  mention  of  further  related  questions  that  deserve  attention 

would  include  (i)  the  (approximate)  distribution  of  max  Q(t),   both  in 

0£t£°° 

sufficient  and  single  server  contexts,  (ii)  approximations  when  the 

basic  process  is  non-Markovian,  (iii)  the  introduction  of  decision 

variables,  e.g.  when  to  begin  servicing,  control  of  arrivals,  etc. 

These  questions,  plus  others,  are  presently  under  active  investigation 

by  the  authors. 
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Table  1:   X=0.25     y  =  1 


N    5    10    15    20     25     30     35     40     45      50 

E[W]   7.6  27.3  68.7  136.4  230.3  350.0  494.6  664.5  859.4  1079.4 
WD    0.1   13.7   54.0   120.6   212.7   330.1  472.5   640.0  832.5   1050.0 


Table  2:   A  =  0.75     U=l 


N     5    10    15     20     25     30     35     40     45      50 

E[W]   10.7  43.3  101.5  184.8  293.1  426.4  584.7   768.1  976.4  1209.7 
W      6.0  36.7   92.5  173.3  279.2  410.0  565.8  746.7  952.5  1183.3 


Numerical  Comparisons  Between  E[W]   and  W 


(Deterministic  Approximation) 
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